function [fl, Pvl, Y] = testspec2(amp, noise_amp, freq, dur)

if (nargin < 1)
   noise_amp = 0;
end

h = findobj('Tag', 'testspec');
if(isempty(h))
   h = figure('Tag', 'testspec', 'name', 'test spectral analysis', 'NumberTitle', 'off');
end;
figure(h);
clf;

fs = 10000; % sample at 10 kHz
RATES = 1/fs; % rate
X=[0:RATES:dur/1000]; % and time in sec
N = randn(size(X));
Y=amp*sin(2*pi*X/(1/freq))+noise_amp*N; % corresponds to 10 msec period = 100 Hz; can add noise too.
subplot('Position', [0.1, 0.55, 0.8, 0.4]);
plot(X,Y);

                  w = hanning(length(Y), 'symmetric');
                  Yp = Y .* w';
%[Pvw, f] = pwelch(Y, 1024, fs, 128, 32);
%[Pvw, f] = pwelch(Y, 1024, fs);
Yf = fft(Yp);
n_pts = floor(length(Y) / 2);
ff = (fs/2) * [0 : n_pts] / n_pts;
Pvf = (abs(Yf(1 : n_pts + 1)) .^ 2) / n_pts;
[Pvw, f] = pwelch(Y, length(X), fs);
Pvw = Pvw(1:end,1)*fs;  % note unscaling by Fs here... to undo what pwelch does...
flomb=[1:300];
[fl, Pvl, conf] = lombpsd2(X, Yp, flomb, 0.95); % compute lomb and get confidence level
Pvw = Pvw/((2 * var(Y))); % rescale fft version in same way as lomb
%Pvf = Pvf/((2 * var(Y))); % rescale fft version in same way as lomb
Pvl = sqrt(Pvl * 2 * var(Yp)/ n_pts) / (sum(w)/length(w));
dur = 1000*(max(X) - min(X));
%Pvl = Pvl / dur;
Pvf = Pvf / dur;
conf = conf / dur;

subplot('Position', [0.1, 0.1, 0.8, 0.4]);
%plot(f(2:end), Pvw(2:end), 'bx-');
%plot(ff, Pvf, 'go-');
hold on;
plot(flomb, Pvl(1:end), 'r-');
plot([flomb(1) flomb(end)], [conf conf], 'k--');
set(gca, 'XScale', 'linear');
set(gca, 'XLim', [min(flomb) max(flomb)]);